TK!2
=v
#1
:n
ne
:v
2
:s
i
:c
number of equations
#2
:n
L
:v
2
:s
i
:u
ft
:d
ft
:c
length of the pendulum
#3
:n
g
:v
32
:s
i
:u
ft/s^2
:d
ft/s^2
:c
gravitational constant
#4
:n
goL
:s
b
:u
1/s^2
:d
1/s^2
#6
:n
w0
:v
0
:s
i
:u
rad/s
:d
rad/s
:c
initial value of angular velocity
#7
:n
th0
:v
.7853981634
:s
i
:u
rad
:d
rad
:c
initial value of angular displacement
#9
:n
x
:v
0
:s
i
:u
s
:d
s
:c
time
:a
x
#10
:n
m
:s
b
:u
rad/s
:d
rad/s
:c
w of equivalent mathematical pendulum 
:a
m
#11
:n
d
:s
b
:u
rad/s
:d
rad/s
:c
diff. vel of math and phys pendulum
:a
d
=u
#1
:f
rad/s
:m
5
:t
rad/s/5
=l
#1
:n
x
:c
independent variable (time)
:v
0
.1
.2
.3
.4
.5
.6
.7
.8
.9
1
1.1
1.2
1.3
1.4
1.5
1.6
E
#2
:n
y
:c
names of dependent variables (functions)
:v
'y1
'y2
'y3
'y4
'y5
'y6
'y7
'y8
'y9
'y10
E
#3
:n
y1
:c
dependent variable (angular displacement)
:v
.7853981634
.729370425099513
.567920165256707
.321713891782996
.0259019639070472
-.273949667105062
-.531357227499456
-.708967001079821
-.783631526339902
-.746462071577142
-.601819437063712
-.367899088206093
-.07751643486847221
.224978901663387
.492488127543551
.6855905428865
.778640108078423
E
#4
:n
y2
:c
dependent variable (angular velocity)
:v
0
-1.10959216887296
-2.08526869011564
-2.78022121745648
-3.0594732355754
-2.85954849565948
-2.22651044603537
-1.28778172491801
-.191228293271535
.927327777452943
1.93549811018881
2.68825920632903
3.04494599698252
2.9258542349393
2.35866318313348
1.46127493504207
.381639888626583
E
#5
:n
k
:c
names of lists with RK coefficients
:v
'k1
'k2
'k3
'k4
'k5
'k6
'k7
'k8
'k9
'k10
E
#6
:n
k1
:c
RK coefficient
:v
1.46127493504207
.954770751978258
.910918613219163
.390320246078474
E
#7
:n
k2
:c
RK coefficient
:v
-10.1300836612762
-11.0071264364581
-10.7095468896359
-11.2146724714648
E
#8
:n
yi
:v
.6855905428865
1.46127493504207
E
#9
:n
ye
:v
.7766824042084171
.390320246078474
E
#10
:n
m
:c
angular velocity, math pendulum
:v
0
-1.22339380337395
-2.25364062518516
-2.9280871453427
-3.14025308804768
-2.85664211605295
-2.12203016579679
-1.05239631157285
.183387780152278
1.39021897372292
2.37756515948554
2.98954608448183
3.129543416895
2.77545465619915
1.98318263200246
.87780967644373
-.366150125741232
E
#11
:n
d
:d
rad/s/5
:u
rad/s
:c
same, difference betw. math and phys pendulum
:v
0
-.113801634500987
-.168371935069515
-.147865927886217
-.0807798524722752
.0029063796065345
.104480280238584
.235385413345165
.374616073423813
.462891196269981
.442067049296726
.301286878152798
.0845974199124751
-.15039957874015
-.375480551131023
-.58346525859834
-.747790014367815
E
=f
#1
:n
RK4
:c
Classical Fourth Order Runge-Kutta method
:t
Procedure
:p
ne
#2
:n
EQ
:c
Differential equations of a physical pendulum
:t
Procedure
:p
goL
:i
t,é,w
:o
dé%dt,dw%dt
#1
:n
RK4
:z
" Notation:  ne   number of 1st order equations
"            EQ   name of the function with 1st order equations
"            x    independent variable (list)
"            y    list with names of dependent variables (lists) y1, y2, ...
"            k    list with names of RK coefficients (lists) k1, k2, ...

" Description: This procedure represents an implementation of a classical
" 4th order Runge-Kutta procedure for numerical integration of differential
" equations.  In this implementation, the set of 1st order equations must be
" available in a procedure function named  EQ .  Whenever this function is
" called from here, the current values of independent variable  x  and
" functions  y1, y2, ... must be passed over as input variables  represented
" by  x , 'ye[1] , 'ye[2] , ... .  The returned values of the right hand
" sides of the 1st order equations are used as RK coefficients  'k1[*] , 
" 'k2[*] , ... (see marked function calls below; they have to be modified
" when the number of 1st order equations changes).

" The values of independent variables must be available in list  x ,  and the 
" solution is accumulated in lists  y1 , y2 , ... .  The names of these lists
" must be declared as symbolic values in the master list  y .  The same
" applies to the lists  k1 , k2 , ... of RK coefficients whose names must
" be declared in the master list  k .  The initial values must be assigned
" as the first elements in the lists  y1 , y2 , ... .

n:= length('x)
for i=2 to n
   xi:= 'x[i-1]
   x:= xi
   for e=1 to ne
       'yi[e]:= 'y[e][i-1]
       'ye[e]:= 'yi[e]
       next e
   h:= 'x[i]-xi
   for j=1 to 3
       call EQ(x,'ye[1],'ye[2];'k1[j],'k2[j])  " <== depends on # of eq-ns
       H:= h/2
       if j=3 then H:=h
       x:=xi+H
       for e=1 to ne
           'ye[e]:= 'yi[e] + H*'k[e][j]
           next e
       next j
   call EQ(x,'ye[1],'ye[2];'k1[4],'k2[4])      " <== depends on # of eq-ns
   for e=1 to ne
       'y[e][i]:= 'yi[e] + dot('k[e],1,2,2,1)*h/6
   next e
next i
call delete('yi)
call delete('ye)
/E
E
#2
:n
EQ
:z
dé%dt:= w                 " definiton of angular velocity
dw%dt:= -goL*sin(é)       " Newton's law applied to physical pendulum

" Note: These are the nonlinear differential equations of a physical pendulum
"       called from the procedure RK4 (classical 4th order Runge-Kutta)
"       which, in turn, is called from the Rule Sheet.
"       For the display of tabulated and plotted solution press  F8  or  F7
"       respectively.  In addition to angular displacement  é  and velocity
"       w  of physical pendulum also displayed are the velocity of equivalent
"       mathematical pendulum, m and five times the difference between the
"       two velocities, d.
/E
E
=r
#1
:r
" *** Differential equations of a physical pendulum - Runge-Kutta method ***
:s
C
#3
:r
goL = g/L
#4
:r
place('y1,1) = th0    " ( initial 
#5
:r
place('y2,1) = w0     " ( conditions
#6
:r
call RK4()
#8
:r
m = -th0*sqrt(goL)*sin(sqrt(goL)*x)   " math pendulum, closed-form solution
#9
:r
d = m - 'y2[elt()]                    " velocity difference, math-phys
:s
C
#11
:r
" Notes:
:s
C
#12
:r
" This model is an example of using the RK4 procedure function for integra-
:s
C
#13
:r
" ting systems of differential equations.  The example deals with equations
:s
C
#14
:r
" of a physical pendulum (see function EQ).  The resulting angular velocity
:s
C
#15
:r
" is compared with the closed-form solution for a mathematical pendulum
:s
C
#16
:r
" obtained by list solving the last two equations above (the  call RK4()
:s
C
#17
:r
" should be cancelled during list solving).  See the Function Subsheets  EQ
:s
C
#18
:r
" and  RK4  for further comments.
:s
C
=p
#1
:n
solution
:t
Line chart
:d
EGA
:i
"Pendulum velocity:  w = phys, m = math, d = 5x difference"
:s
Yes
:z
X-axis
:g
No
:l
None
:h
"x"
:v
""
:x
x
:u
#1
:y
y1
:f
1
:p
Both
:c
"é"
#2
:y
y2
:f
1
:p
Both
:c
"w"
#3
:y
m
:f
1
:p
Both
:c
"m"
#4
:y
d
:f
1
:p
Both
:c
"d"
E
=t
#1
:n
pendulum
:s
Screen
:v
Vertical
:f
1
:r
""
:c
""
:t
"Physical Pendulum"
:u
#1
:l
y1
:w
10
:h
""
#3
:l
y2
:w
10
:h
""
#5
:l
m
:w
10
:h
""
#6
:w
10
:h
""
#7
:l
d
:w
10
:h
""
E
=n
=g
#1
:c
.000001
:i
10
:v
Yes
:l
60
:f
66
:m
0
:w
80
:b
Yes
:a
Yes
:r
Yes
:h
No
:t
1
:n
Yes
:d
LPT1
%Tf,f,0,0,0,1,0>t,0,1,0,5,0
